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ABSTRACT 

We reexamine both numerically and analytically the collapse of the singular isothermal sphere in the con- 
text of low-mass star formation. We consider the case where the onset of collapse is initiated by some arbi- 
trary process which is accompanied by a central output of either heat or kinetic energy. We find two classes of 
numerical solutions describing this manner of collapse. The first approaches in time the expansion wave solu- 
tion of Shu, while the second class is characterized by an ever-decreasing central accretion rate and the pres- 
ence of an outwardly propagating weak shock. The collapse solution which represents the dividing case 
between these two classes is determined analytically by a similarity analysis. This solution shares with the 
expansion wave solution the properties that the gas remains stationary with an r~ 1 2 density profile at large 
radius and that, at small radius, the gas free-falls onto a nascent core at a constant rate which depends only 
on the isothermal sound speed. This accretion rate is a factor of ~0.1 that predicted by the expansion wave 
solution. This reduction is due in part to the presence of a weak shock which propagates outward at 1.26 
times the sound speed. Gas in the postshock region first moves out subsonically but is then decelerated and 
begins to collapse. The existence of two classes of numerical collapse solutions is explained in terms of the 
instability to radial perturbations of the analytic solution. Collapse occurring in the manner described by 
some of our solutions would eventually unbind a finite-sized core. However, this does not constitute a vio- 
lation of the instability properties of the singular isothermal sphere which is unstable both to collapse and to 
expansion. To emphasize this, we consider a purely expanding solution for isothermal spheres. This solution is 
found to be self-similar and results in a uniform density core in the central regions of the gas. Our solutions 
may be relevant to the “ luminosity ” problem of protostellar cores since the predicted central accretion rates 
are significantly reduced relative to that of the expansion wave solution. Furthermore, our calculations indi- 
cate that star-forming cloud cores are not very tightly bound and that modest disturbances can easily result in 
both termination of infall and dispersal of unaccreted material. 

Subject headings: hydrodynamics — shock waves — stars: formation 


1. INTRODUCTION 

The gravitational collapse of gas spheres in the context of 
star formation has long been an active area of study. Pioneer- 
ing numerical calculations (e.g., Gaustad 1963; Bodenheimer 
& Sweigart 1968; Larson 1969; Penston 1969a, b) based on the 
collapse of a Jeans unstable uniform density sphere elucidated 
the general properties of star formation and provided the basis 
of much of the subsequent work in the area. One of the most 
important properties found by these early investigations was 
that efficient cooling by dust grains allowed the rapid radiation 
of compressional heat generated during collapse. The gas 
remained essentially isothermal with a temperature of ~ 10 K 
over many orders of magnitude in density during the early 
phases of collapse. Isothermality was only violated when the 
densities became sufficiently high for gas to become rather 
optically thick to dust emission. Isothermality, taken as an 
assumption, enabled great simplification in many subsequent 
studies because radiation transfer was eliminated from the 
problem (however, this restricted the applicability of these 
solutions to the more optically thin regions of collapse). It was 
also shown that collapse occurred in a very nonhomologous 
fashion, i.e., the central regions of the cloud collapsed much 
more rapidly than the outer regions. The density distribution 


inevitably approached a strongly peaked r ~ 2 profile as the 
collapse approached core formation. 

Later studies addressed primarily the unrealistic nature of 
the starting conditions adopted in the early work. The assumed 
initial configuration for the gas was an isothermal uniform 
density sphere that was far out of hydrostatic equilibrium. 
Such conditions would be very difficult to produce in nature. 
Furthermore, the velocities resulting from the collapse were 
criticized as being unrealistic (Shu 1977). At large radius, the 
flow was directed inward at 3.3 times the sound speed which 
was considered numerically ad hoc. Many of the subsequent 
considerations (e.g., Shu 1977; Hunter 1977; Boss & Black 
1982; Whitworth & Summers 1985; Foster & Chevalier 1993) 
centered around the collapse of equilibrium isothermal spheres 
(Bonner-Ebert spheres; see Ebert 1955; Bonner 1956) where 
the gas was assumed to be isothermal at all times and initially 
stationary. Because the initial configurations were in equi- 
librium, the numerical calculations of collapse were initiated 
by some rather ad hoc perturbation. However, the nature of 
the collapse was found to be rather insensitive to the manner of 
the initial perturbation (e.g., Hunter 1977). 

Many collapse solutions were found. These differed from one 
another based on the initially assumed density profile and 
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ranged from the collapse of the marginally unstable Bonner- 
Ebert sphere (Foster & Chevalier 1993) to the maximally 
unstable singular isothermal sphere (Shu 1977; Hunter 1977). 
However, it was not clear which, if any, was the most relevant 
to present-day low-mass star formation. Recently, both theo- 
retical and observational work has yielded a paradigm for low- 
mass star formation (see Shu et al. 1993 and references therein). 
A cloud initially supported against collapse by gas pressure, 
turbulence, and magnetic fields evolves quasistatically toward 
a centrally condensed singular isothermal sphere due to the 
outward diffusion of the magnetic field (ambipolar diffusion). 
After ^10 6 yr, the distribution becomes sufficiently centrally 
condensed to undergo an “inside-out” collapse in the manner 
of the analytic, self-similar solution of Shu (1977) (dubbed the 
“expansion wave solution,” hereafter EWS). Although the col- 
lapse calculation has been generalized to account for modest 
amounts of rotation (Terebey, Shu, & Cassen 1984) and mag- 
netic fields (Galli & Shu 1993a, b), the simplicity of the original 
solution, as well as its closeness to the generalized solutions at 
relatively early times and in regions somewhat outside the 
core, has led it to be generally adopted as the benchmark for 
testing the agreement of the paradigm with observations. 

The EWS is obtained when a static singular isothermal 
sphere suffers a central perturbation which causes the core to 
collapse. The ensuing collapse proceeds in an “inside-out” 
fashion where the infalling region is bounded by a rarefaction 
wave which propagates outward at the sound speed. Outside 
the rarefaction wave, the gas remains undisturbed and in 
hydrostatic equilibrium. At small radius, gas free-falls onto a 
growing core at a constant rate which depends only on the 
sound speed a and the gravitational constant G, 

M = C~, (1) 

where C = 0.975 for the EWS. The collapse can be initiated by 
any process which causes an unbalancing of the two forces in 
the core (Shu 1977) or by the onset of true dynamical insta- 
bility (Shu, Adams, & Lizano 1987). In the former case, any- 
thing which decreases the central gas pressure, increases the 
gravity (mass), or both, will cause core collapse. Possible physi- 
cal processes include molecule formation (Shu 1977) and 
primary core collapse (e.g., Larson 1969; Winkler & Newman 
1980). True dynamical instability arises during the evolution of 
the protostellar cloud toward ever more centrally peaked 
density profiles during the process of ambipolar diffusion. 

The considerations of this paper derive from the observation 
that the suggested processes which might initiate the EWS 
collapse will not only cause an initial unbalancing of the 
opposing central forces, but will also in general be accompa- 
nied by some initial or sustained burst of energy which will be 
reabsorbed by parts of the collapsing cloud. Certainly, imme- 
diately after the onset of collapse in the inside-out picture, 
energy will be released from a radiating accretion shock 
around either a core or a disk, or possibly by gas streams, 
which are not flowing exactly radially, colliding outside of the 
core. This energy will propagate back into static and sub- 
sonically moving parts of the envelope and cause a preheating 
effect. That is, the radiation will heat the gas above the 
assumed isothermal temperature, causing an outward push 
due to a higher gas pressure. This extra pressure may have 
significant dynamical effects. Other ways in which collapse can 
be initiated may also be accompanied by energetic output, 


such as molecule formation which produces heat via recombi- 
nation. Additionally, sustained stellar winds produced in the 
very early life of a nascent protostar could push material out 
via kinetic pressure. In any case, since the conditions which 
initiate protostellar collapse are poorly understood, it is impor- 
tant to investigate a spectrum of possibilities, the results of 
which may be compared to observations in order to better 
understand star formation. 

We can imagine intuitively what the effects might be of 
absorption of an initial burst of energy followed by an 
unbalancing of the pressure in the core of a static, isothermal 
sphere. The extra energy or pressure will first cause an outward 
push to the gas in the central regions. In one instance, if the 
amount of energy or the outward push is very small, the initial 
condition may be thought of as a perturbation on the condi- 
tion which gives rise to the EWS. We might imagine that in 
time, as the ensuing collapse progresses to encompass scales 
much larger than that set by the initial condition, the effects of 
the perturbation will die out, and the EWS will be recovered. 
In another case, if the outward push is stronger, a shock will 
likely form. Because of the steep density gradient of the singu- 
lar isothermal sphere, this shock can coast to large radius. In 
the center, after the initial burst has subsided, gas which was 
pushed outward will decelerate due to the pressure imbalance, 
eventually going into infall. 

We numerically determine collapse solutions resulting from 
the energy burst initial condition in § 2. However, in the case 
where a shock forms, we might expect in analogy with the 
Sedov-Taylor solution (e.g., Sedov 1959) that there would be 
an evolution toward a similarity solution once the shock has 
moved to distances much greater than the scale set by the 
initial perturbation. Recall that the Sedov-Taylor solution 
describes the state of gas, for example, after a piston pushes 
into an initially static uniform density gas (gravity is not 
considered). The motion of the piston gives rise to a shock. 
When the shock has propagated into a region far removed 
from the scale set by the motion of the piston, the solution 
approaches self-similarity because the initial length scale 
becomes irrelevant. By analogy, when the shock in our 
problem has moved beyond the core area where the initial 
energy burst occurred, there are also no longer any relevant 
length scales. Contrary to this view, we find that the numeri- 
cally generated solutions do not tend toward any similiarity 
solution; rather, they depend on the particular conditions initi- 
ating the collapse. However, this result can be understood in 
terms of the instability of a self-similar solution which we 
present in § 3. In § 4, we introduce an expanding solution to 
help elucidate the properties of the new collapse solutions as 
well as the general collapse properties of the singular isother- 
mal sphere. We discuss the relation of the new solutions to the 
EWS and the theory of star formation in § 5. Conclusions are 
also given in § 5. 

2. NUMERICAL SOLUTION 
2.1. Method 

The case where radiation from an accretion front preheats 
gas which is not yet in a state of free fall is best treated by a 
hydrodynamical code with radiation transfer. Although similar 
calculations of collapse have been carried out (e.g., Morfill, 
Tscharnuter, & Volk 1985; Bodenheimer et al. 1990; Yorke, 
Bodenheimer, & Laughlin 1993), the initial conditions were 
always assumed to be somewhat out of hydrostatic equi- 
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librium. We will report elsewhere on a one-dimensional calcu- 
lation treating the effects of accretion luminosity assuming the 
exact singular isothermal sphere as the initial state. In the 
present paper, we take a less specific approach by approx- 
imating the effects of the energy deposition by a piston. 
Although this approach is more heuristic than solving a radi- 
ative transfer problem, it lends itself to a general elucidation of 
the relevant physics of the problem and the consideration of 
other possible modes of energy deposition in the center of the 
cloud. 

We assume that the protostellar cloud is spherically sym- 
metric, is always isothermal, and obeys the ideal gas equation 
of state. The dynamical calculations are carried out with a 
Lagrangian hydrodynamics code modified from that devel- 
oped by Yahil, Johnston, & Burrows (1987). The code is con- 
servative in the fluxes of energy and momentum and has 
second-order spatial accuracy in the fluxes. 

The initial configuration of the cloud is a singular isothermal 
sphere. We begin a calculation by assuming that the material 
within a radius of 2 x 10 13 cm, which has an enclosed mass of 
2.5 x 10 -4 M 0 for a cloud temperature of 20 K and a mean 
molecular weight of 2, forms a protostellar core. If there is no 
energy input, and the material immediately outside of this 
region is forced into free fall, we recover the EWS. 

A major problem with using a Lagrangian code for an accre- 
tion calculation is that the Courant time can be so reduced for 
the region near the central object that it makes the investiga- 
tion of the long-term evolution of the rest of the object 
extremely difficult. We circumvent this problem by specifying a 
critical radius such that if a shell of material (corresponding to 
a zone in the calculation) falls below this critical radius, it is 
assumed to have been accreted onto the protostellar core, 
which then only acts gravitationally on the rest of the cloud. 
The mass accretion rate can then be approximated by m shelJ /Af, 
where m shell is the mass of the shell and At is the time elapsed 
since the last shell was accreted. This accretion rate is then held 
constant until the next shell is accreted. The accretion rate 
calculated this way for the EWS follows the analytical value 
closely, with < 1% variations, due to changes of zone sizes and 
time steps. 


We simulate the input of the initial burst of energy by apply- 
ing a piston with a velocity v p immediately outside the proto- 
stellar core and starting at the critical radius of 2 x 10 13 cm. 
The piston is pushed at a constant velocity until it has reached 
some distance d p , where the piston is instantaneously stopped 
and we apply the free-fall boundary condition. This boundary 
condition corresponds to the limiting case where the inner 
pressure is reduced to zero in order to initiate collapse. Tests 
with different boundary conditions where the piston is allowed 
to decelerate gradually under the influence of less severe pres- 
sure gradients give similar results, so they are not discussed. 
We assume a constant pressure outer boundary condition. 

2.2. Results 

In the first case considered (case [i]), we move the piston at 
half the isothermal sound speed (a = 2.87 x 10 4 cm s" 1 for our 
chosen conditions) to a distance of d p = 2 x 10 14 cm. The 
resulting velocity and density profiles of the gas at different 
times are plotted in Figures la and 2a, respectively, and the 
accretion rate at the core radius is given by the uppermost 
solid curve of Figure 3. The action of the piston forms a shock 
which propagates outward. Despite the subsonic motion of the 
piston, a shock is formed because the buildup of material in 
front of the piston causes a supersonically moving front to 
advance ahead of the piston. The material outside of the shock 
remains stationary with a density given by that of the singular 
isothermal sphere. Inside the shock, gas is flowing outward 
subsonically, having been accelerated by the passage of the 
shock. The gas eventually is decelerated by the pressure deficit 
at small radius and goes into free fall at the inner boundary. 
The strength of the shock declines with time and eventually 
becomes an acoustic wave propagating outward at the sound 
speed. The region of outflowing gas and the outward velocity 
of postshock gas become correspondingly smaller as the shock 
degrades. Similarly, the size of the discontinuity in the density 
profile decreases with time. The accretion rate starts out low 
(after settling down from transients) but after a very short time 
rises rapidly and tends toward the EWS value appropriate to 
the chosen temperature of 20 K. 

This case corresponds to a collapse which is initially per- 
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Fig. 1 a Fig. 16 

Fig. 1.— (a) The velocity profiles of the gas for the case where v p - 0.5a and d p = 2 x 10 14 cm are plotted at several different times. The leftmost curve 
corresponds to the earliest time of 0.1 x 10 12 s after starting the collapse. Moving toward the right, the curves correspond to times of (in units of 10* 2 s), 0.307, 0 507, 
0.723, 0.967, 1.25, 1.57, 1.97, 2.42, 2.95, 3.58, 4.35, and 5.24. (6) The velocity profiles of the gas for the case where v p = 0.5a and d p = 2 x 10' 5 cm are plotted at times 
(in units of 10 12 s) 0.239, 1.18, 2.63, 4.74, 7.6, 1 1.6, and 16.7. The leftmost curve gives the earliest time, and the rightmost curve gives the latest time. 
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Fig. l.—ia) The density profiles of the gas for the case where v p - 0.5a and d p - 2 x 10 14 cm are plotted for the same set of times as Fig. la. The topmost curve 
corresponds to the earliest time of 0.1 x 10 12 s, and successively lower curves correspond to later times, (b) The density profiles of the gas for the case where = 0.5a 
andd, = 2 x 10 15 cm are plotted for the same set of times as Fig. 1 b. The topmost curve gives the earliest time, and lower curves give later times. 


turbed from the EWS by the action of the piston. (Recall that 
we obtain the EWS by simply unbalancing the forces in the 
center with no accompanying piston moving outward.) 
However, the effects of the perturbation die out rapidly, and 
the collapse approaches the EWS solution. 

On pushing the piston out farther to d p = 2 x 10 15 cm with 
the same velocity v p (case [ii]), we obtain the gas velocity and 
density profiles shown in Figures 1 b and 2b, respectively, and 
an accretion rate given by the lowest solid curve of Figure 3. 
We again get a shock formed by the early action of the piston, 
but the shock persists in time with constant amplitude all the 
way to the outer boundary. As before, the gas is stationary 
outside the shock and outflowing immediately inside of it. 
However, because of the undiminished amplitude of the shock, 
the outflowing gas occupies a significant fraction of the 
nonstatic region. The shock is rather weak with a Mach 
number of - 1.26, and the gas has a postshock outflow velocity 



Fig. 3. — The central mass accretion rates for several cases are plotted vs. 
time. All cases shown ( solid lines ) assume v p = 0.5a but have different values of 
d p . Starting from the topmost curve and going toward the lowest curve, the 
assumed values of d p are 2 x 10 14 cm, 4 x 10 14 cm, 6 x 10 14 cm, 7 x 10 14 cm, 
and 2 x 10 15 cm. The upper dotted line gives the constant mass accretion rate 
of the EWS (see eq. [1]), and the lower dotted line gives the constant rate of the 
shock solution. 


of ^ 0.43 a. The postshock density rises above the value given 
by the singular isothermal sphere for a considerable region 
inside the shock. The central mass accretion rate starts out 
very low and declines with time and never rises back up to the 
EWS value. In pushing the piston out to larger distances than 
in case (i), we have essentially perturbed the collapse farther 
away from the EWS solution. We might have expected that the 
net effect would be for the collapse to take a much longer time 
than that of case (i) to evolve back to the EWS. Instead, the 
collapse does not show a tendency to evolve back to the EWS 
at all. 

Several cases with values of d p between those of cases (i) and 
(ii) above are also computed (v p = 0.5a in all these cases). The 
central accretion rates for these collapses are shown in Figure 
3. The cases where the central accretion rate increases with 
time have gas velocity and density profiles which behave in 
time similar to those of case (i), and cases with declining accre- 
tion rates have velocities and densities similar to those of case 
(ii). We see that two kinds of solutions are obtained from the 
piston experiments. When the collapse is initiated by a piston 
which pushes to something less than ~6 x 10 14 cm, the nature 
of the subsequent flow is to evolve secularly back toward the 
EWS. For simplicity we call these class I solutions. Once the 
piston is allowed to push farther out than this, the evolution of 
the flow is toward ever decreasing accretion rates (class II 
solutions). The farther out the piston is allowed to go, the 
smaller is the accretion rate. 

In Figure 4, we show the mass of the inner sink cell (the 
“ protostar ”) as a function of time for the cases shown in 
Figure 3. The times required to accrete a sizable mass ( ~ 1 M 0 ) 
onto the sink cell can be significantly extended relative to the 
dynamical time if the collapse proceeds according to the lower 
curves of Figure 4. In these cases, the mass of the sink cell 
would only be ~0.1 M Q by the time the shock hits the bound- 
ary of the calculation. If star formation were to occur in this 
manner, the accretion of a large fraction of the mass of the 
protostar occurs after the passage of the shock through the 
boundary, assuming a typical size for the cloud core. The sub- 
sequent collapse would then depend on the particular proper- 
ties of this boundary. 

We obtain the same set of results as above for different 
assumed piston speeds. We do not show the gas velocities and 
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Fig. 4. — The mass of the central star is plotted as a function of time for the 
same cases as shown in Fig. 3. The top curve corresponds to the case where 
d p = 2 x 10 14 cm, the next highest curve gives the case where d p = 4 x 10 u 
cm, and so forth for the lower curves. The dotted line gives the mass for the 
EWS solution. 


densities in detail. Figure 5 shows the mass accretion rates for 
some of the cases considered. For various piston speeds, the set 
of collapse solutions again divides into those cases with rising 
accretion rates and those with declining rates. The dividing line 
between these two sets of solutions occurs for d p ~ 6 x 10 14 
cm for the piston velocities considered. 

We examine the above results in more detail. The condition 
which determines which of the two classes of solutions a given 
collapse calculation will follow is not easily defined. First, the 
amount of work done by the piston during the push phase of 
the calculation will not alone determine the subsequent evolu- 
tion of the collapse. Having taken the gas to be isothermal 
implicitly assumes that the cloud can instantaneously 
exchange energy with some external sources and sinks. As the 
piston pushes the gas outward, compressed regions of the gas 
will radiate away energy. Thus, the change in the energy of the 
cloud during the time that the piston is pushing does not 
depend solely on the amount of work done by the piston, but 



also on its exact velocity history. Also, the collapse after the 
piston has been stopped is not energy conserving, again since 
the gas is assumed to be isothermal. Whatever changes in the 
energy of the cloud occurred during the time the piston was 
pushing are swamped by subsequent energy exchanges of the 
cloud with the external sources and sinks. 

An alternative is to consider the dynamics of the shock. The 
shock is launched by the piston and grows in strength as the 
piston pushes outward. However, when the piston is stopped 
and the inner boundary condition is applied, a rarefaction 
wave is also sent propagating into the gas because the velocity 
of the gas just in front of the piston is positive when the piston 
is stopped. This can be seen in the topmost curve 
(corresponding to the earliest time shown) of Figure 2b. The 
shock is at a radius of ~10 16 cm, and a rarefaction wave is 
seen at a radius of ~6.3 x 10 15 cm. A rarefaction wave travels 
at the local sound speed. Since the rarefaction wave generated 
by the piston propagates into a region of outflowing gas, the 
wave can actually catch up to the shock, even though the shock 
is traveling outward supersonically. This can be seen in both 
Figures 1 b and 2b. The rarefaction wave which is well separat- 
ed from the shock at the earliest time plotted has caught up to 
the shock by the next time shown and is interacting with the 
shock and decreasing its strength. Once caught up with the 
shock, the rarefaction wave cannot proceed further upstream 
since the shock is propagating supersonically into a static 
medium. The subsequent evolution of the shock and the col- 
lapse solution in general then depends on the complicated non- 
linear interactions of the rarefaction wave and the shock. The 
amplitudes of the shock and the rarefaction wave again depend 
in detail on the velocity history of the piston. 

Although complex, there is hope that the late time behavior 
of the collapse can be understood. The evolution toward the 
EWS of collapse initiated by central perturbations without 
accompanying piston motions suggests that there must be 
some link between the late time, scale-free evolution of our 
collapse solutions with a similarity solution. Analogous con- 
siderations based on the Sedov-Taylor problem further 
emphasize this link. What then should be the properties of 
such a similarity solution? The numerical considerations indi- 
cate that if the actions of the piston are very finely tuned, a 
collapse which has properties between those of the two 
numerical classes can be obtained. Specifically, the solution 
would have a constant central accretion rate and an approx- 
imately Mach 1.26 constant amplitude shock propagating 
outward. We find such a solution in § 3 and use it to under- 
stand the late time evolution of the numerical solutions. 


3. SIMILARITY SOLUTIONS 
3.1. The Shock Solution 

The equations governing an isothermal gas sphere are the 
continuity equation. 


dp 

dt 


II 

r 2 dr 


( r 2 pu ) = 0 , 


( 2 ) 


and the Euler equation. 


Fig. 5. — The central mass accretion rates for several cases are plotted as in 
Fig. 3. The solid lines give the case where v p - a. The upper solid line assumes 
d p = 4.7 x 10* 4 cm, and the lower one assumes d p — lA x 10 14 cm. The 
dashed curves give the case where u = 1.2a. The upper curve assumes d p = 5.9 
x 10 14 cm, and the lower curve assumes 6.3 x 10 14 cm. The dotted lines give 
the accretion rates of the EWS (top line) and the shock solution (bottom line). 


du du 

at dr 


dp 

dr 


GM 

„2 


(3) 


where p is the gas density, u is the gas velocity, M is the mass 
enclosed within radius r, and a is again the isothermal sound 





No. 2, 1995 


PROTOSTELLAR COLLAPSE WITH SHOCK 


779 


speed. Following Shu (1977), we define the similarity variable 


r 

x = — , 
at 

and we nondimensionalize the fluid variables, 


(4) 


P(r, t) = 


«(*) 
4nGt 2 ’ 


A/(r, t) = — m(x ) , 

o 


(5) 

( 6 ) 


u(r, t) = ar(x) . 


(7) 


Equation (6) indicates that the mass is zero at the instant t = 0. 
This corresponds to the moment of core formation. We will 
consider positive times after the formation of the core. Substi- 
tuting equations (4)-(7) into equations (2) and (3), we get the 
following set of equations (same as eqs. [11] and [12] of Shu 
1977), 


[(x - v) 2 



a(x — r) 


2 

x 


(X - V) , 


( 8 ) 


[(x - v) 2 - 1] 


1 da 
a dx 


a — - (x — i?) 
x 


(x - v) , 


(9) 


for the nondimensionalized fluid variables. It can also be 
shown using the continuity equation for M(r, t) rather than for 
the density that the nondimensionalized mass is given by 

m = x 2 a(x — i>) . (10) 


Shu (1977) found a large set of solutions to equations (8)- 
(10), some of which were identified as viable collapse solutions. 
These were called “minus solutions without critical points” 
(MSWCPs) (see Fig. 2 of his paper and the light dotted lines of 
Figs. 6 and 7), of which the EWS represents a certain limit 
(although the EWS is not one of the MSWCPs). The MSWCPs 
correspond to solutions of equations (8)-(10), where the outer 
boundary condition requires that the gas be static at large x 
and the inner boundary condition requires that the gas be in 
freefall, 

m -*> m 0 , a -> (m 0 /2x 3 ) l/2 , v — (2 m 0 /x) 1/2 as x -► 0 , 

( 11 ) 


where m 0 is the nondimensionalized mass of the accreting 
central protostar. These solutions are so dubbed because the 
entire flow from the outer static region to the inner free-falling 
regime occurs without passing through a critical point. Note 
that the coupled set of equations (8) and (9) have a locus of 
critical points given by 

x=l+», (12) 

which is shown as the dashed line of Figure 6. The entire set of 
MSWCPs does not intersect this line. 

The MSWCPs correspond physically to the collapse of ini- 
tially static isothermal spheres with r 2 density profiles but 
which are too dense to be in equilibrium. The initial density 
can be represented by 

P(r,0) = ^r- 2 , (13) 


where A is a numerical constant satisfying the condition A > 2. 



X 

Fig. 6. — The nondimensionalized velocities for various similarity solutions 
are plotted vs. the similarity variable. The light dotted lines give the “minus 
solutions without critical points” (MSWCPs). From the lowest to the upper- 
most curve, these correspond to A = 2.8, 2.6, 2.4, and 2.2. Hie heavy dotted 
line gives the expansion wave collapse solution (EWS; A — 2), and the heavy 
solid line gives the shock solution. The light solid line gives the expanding 
solution, and the dashed line gives the locus of critical points of eqs. (8) and (9) 
of the text. 



-2 -1.5 -1 -.5 0 .5 

log X 


Fig. 7. — The nondimensionalized density profiles for several similarity 
solutions are plotted. The light dashed line gives the MSWCPs corresponding 
to A = 2.8, and the heavy dotted line gives the expansion wave collapse solu- 
tion (EWS; A = 2). The heavy solid line gives the shock solution, and the light 
solid line gives the expanding solution. 
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The limiting case of A = 2 corresponds to the singular isother- 
mal sphere, implying that the EWS is the limit of the MSWCP 
solutions as A -> 2. The values of m 0 are determined uniquely 
by the value of A. For the EWS, m 0 = 0.975, and it is larger for 
solutions with larger values of A. 

The EWS obeys the same inner and outer boundary condi- 
tions as the MSWCPs; however, the flow passes through a 
critical point at x = 1 (Fig. 6). This implies that the EWS is 
actually a member of a different set of solutions called “minus 
solutions” (Shu 1977). These are solutions which obey the 
inner boundary condition of equation (11) (with values of 
m 0 < 0.975) but which do pass through a critical point with 
positive slope in the velocity profile. With the exception of the 
EWS, the minus solutions were originally discarded as inap- 
propriate for collapse since the inner inflowing part of the 
solution could not be joined consistently with an outer solu- 
tion which obeyed the static outer boundary condition. More 
specifically, starting with values given by the inner boundary 
condition (1 1), the fluid variables can be integrated outward to 
the location of the critical line (e.g., the negative velocity 
portion of the heavy solid line in Fig. 6 is a minus solution). In 
order that the flow then passes smoothly through the critical 
point with positive slope in the velocity profile, equations (8) 
and (9) require that 


-U = (1 - x ) (x-x*) + ••• 

** 

2 2 

a = — - -j (x - x J + • • ■ , 


(14) 


where x * is the position of the critical point (the value of x m is 
fixed by the value of m 0 or vice versa). The fluid variables can 
then be integrated outward to the point where the gas velocity 
is zero. Here the velocity profile of the inner solution can be 
smoothly joined to the static outer part. However, the density 
at this point from the inner integration does not correspond to 
the density appropriate to that radius in the singular isother- 
mal sphere. The velocity and density therefore cannot be 
simultaneously matched to a static outer solution. The EWS is 
the only exception that allows a consistent matching. 

Guided by the numerical considerations of § 2, the presence 
of a shock somewhere in the flow allows for 1 more degree of 
freedom (the location of the shock) which might allow the fluid 
variables of the inner part of a minus solution to be consistent- 
ly joined via the shock jump conditions to a static outer solu- 
tion. We first determine the appropriate shock jump 
conditions. Since we consider an isothermal shock, we need 
only consider mass and momentum conservation across the 
shock. If we assume that the shock is propagating into static 
gas, we have, in terms of dimensional fluid variables, 





(15) 


where u 2 and p 2 are the postshock gas velocity and density, 
respectively, u 9h is the speed of the shock as it goes into the 
static gas, and p, is the preshock gas density. Assuming the 
shock travels at a constant speed and starts at the origin at 
t = 0, the location of the shock expressed in terms of the simi- 
larity variable will be fixed. That is, u sh = x sh a , where x sh is the 
fixed position of the shock in similarity space. The jump condi- 
tions in nondimensional form are then 



(16) 


where again a subscript of 2 denotes postshock values. In 
deriving equation (16), we have taken the preshock density to 
be that of the singular isothermal sphere. 

The property that a 2 = 2 at the shock suggests the following 
scheme for determining the desired shock solution. Pick a loca- 
tion on the line of critical points such that x* lies in the range 
0 < x* < 1. The values of the fluid variables at x* are then 
fixed by equation (14). Using these as starting values for an 
integration to smaller x yields the value of m 0 which corre- 
sponds to the assumed value of x*. The fluid variables are then 
integrated outward from x* to the point where a 2 = 2. Take 
this as the location of the shock x sh . The jump condition for the 
density is then automatically satisfied at this point. However, 
the velocity jump condition may not be simultaneously satis- 
fied. The initially assumed value of x* is adjusted and the 
above procedure repeated until both jump conditions are 
simultaneously satisfied. 

We find that there is a unique value of x* (hence also of m 0 
and x sh ) such that the solution both passes smoothly through 
the critical point and satisfies the shock jump conditions, given 
that the outer solution is the singular isothermal sphere. This 
solution is shown as the heavy solid line in Figures 6 and 7, and 
a listing of values is given in Table 1. We get x* = 0.0544, 
which implies that x sh = 1.26 and that m 0 = 0.105. That is, the 
shock has a Mach number of 1.26, and the mass accretion rate 
onto the central “protostar” is constant and is given by equa- 
tion (1) with C = 0.105. 

This solution represents the dividing line between the 
numerical class I and class II solutions. The predicted constant 
accretion rate neatly separates those cases with increasing acc- 


TABLE I 


Shock Solution* 


X 

( 1 ) 

a 

( 2 ) 

V 

( 3 ) 

m 

( 4 ) 

0.00 

00 

— QO 

0.105 

0.05 

40.0 

- 1.03 

0.108 

0.10 

21.2 

- 0.446 

0.116 

0.15 

15.7 

- 0.218 

0.130 

0.20 

13.0 

- 0.0921 

0.152 

0.25 

11.3 

- 0.00891 

0.182 

0.30 

10.0 

0.0532 

0.222 

0.35 

9.02 

0.103 

0.272 

0.40 

8.19 

0.146 

0.333 

0.45 

7.47 

0.183 

0.403 

0.50 

6.83 

0.217 

0.484 

0.55 

6.26 

0.247 

0.574 

0.60 

5.75 

0.275 

0.673 

0.65 

5.29 

0.300 

0.781 

0.70 

4.86 

0.324 

0.896 

0.75 

4.48 

0.345 

1.02 

0.80 

4.13 

0.365 

1.15 

0.85 

3.80 

0.383 

1.28 

0.90 

3.51 

0.400 

1.42 

0.95 

3.24 

0.414 

1.57 

1.00 

2.99 

0.427 

1.71 

1.05 

2.77 

0.439 

1.87 

1.10 

2.56 

0.448 

2.02 

1.15 

2.37 

0.455 

2.18 

1.20 

2.19 

0.461 

2.33 

1.25 

2.03 

0.463 

2.49 


* Col. (!) gives the position in similarity 
space, col. (2) is the nondimensionalized 
density, col. (3) is the nondimensionalized 
velocity, and col. (4) is the nondimension- 
alized enclosed mass. 
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retion rates from those with declining rates (see Figs. 3 and 5). 
In addition, since the analytic “shock ” solution is self-similar, 
the distance of the shock to the origin is a fixed multiple of the 
distance of the stagnation point to the origin. This property 
separates the class I solutions which have a decreasing value of 
the ratio of these two distances (Fig. la) from the class II 
solutions which have an increasing ratio (Fig. lb). The charac- 
ter of the shock in the analytic solution is also comparable to 
that of the class II solutions, moving out at about the same rate 
and having constant amplitude. 


3.2. The Energetics of the Shock Solution 

We consider in more detail some of the properties of the 
shock solution. An obvious issue is what drives the large region 
of outflowing gas. This is most easily addressed by looking at 
how the energy of the solution evolves with time. In principle, 
the static singular isothermal sphere has infinitely negative 
total energy. This is because the total energy is the sum of the 
gravitational energy (E grav ) and the thermal energy (£ th ), which 
are each given within a radius r by 


£.r.vM = “ 


p(r)47ir 2 dr 


j r GM( 



(17) 

(18) 


where gm p is the mean mass per particle. Clearly, if r extends to 
infinity, the total energy is infinitely negative. Note that equa- 
tions (17) and (18) are consistent with the virial theorem, since 
the singular isothermal sphere truncated at any finite radius is 
in equilibrium only if an external pressure acts at the outer 
boundary. A boundary term, which is customarily discarded in 
the statement of the theorem, must be included in our case. 
Consider now the energy of the flowing part of the shock solu- 
tion. Within a fixed value of the dimensionless distance x, 


^grav(^’ 0 


a 5 t 

~G 


*x 

xct(x)m(x)dx , 

Jo 


3 a 5 t _ . . - 

£,„(*, 0 = 2 ~c X m °l ’ 


(19) 

( 20 ) 


where the fluid variables are converted to the forms of equa- 
tions (5) and (6). The quantity m 0 in equation (20) represents 
the mass accreted onto the core and is subtracted from the 
total enclosed mass, since the accreted gas is assumed to have 
no thermal effect on the collapse. We additionally need to 
consider the total kinetic energy of the flowing gas, 

£ ki „(x) = | cc(x)u 2 (x)x 2 dx . (21) 

The integrals of equations (19) and (21) can be evaluated 
numerically within the shock radius x sh . The total energy is 
then 


^tot(*sh> 0 — ^grav F E t h 


(-3.53 + 3.62 + 0.17)^ = 0.26^, (22) 
G G 


which is positive. This implies that as the flowing region 
expands outward in time, the energy of each region that gets 
engulfed by the expanding shock increases. This behavior can 
be contrasted with that of the collapsing region of the EWS. 


Equations (19)— (21) imply that inside of the rarefaction wave 
x = 1, the energy is 

EJLx = 1, t) Egrav "1” ”1” £kin 

= ( — 2.98 + 1.54 + 0.69) ^ = —0.75 -p - . (23) 
G G 

We consider the flow of energy in these solutions. The iso- 
thermal nature of the gas does not derive from the gas being in 
thermodynamic equilibrium with a large heat reservoir. In the 
realistic setting of modeling the collapse of a cloud core, iso- 
thermality arises out of the stable balance between heating due 
to sources such as cosmic rays and cooling due to radiation 
from dust. For example, in the cases where the gas is heated by 
compression, dust emission is assumed to be so efficient that 
the extra heat can be radiated away rapidly. In the opposite 
case where gas is cooled by expansion, the various sources of 
heating are assumed to be sufficiently rapid to heat the gas 
back to the stable equilibrium temperature. In this sense, the 
shock solution does not represent a violation of the second law 
of thermodynamics. That is, if we had considered a cloud 
which was isothermal with temperature T because it was in 
thermal contact with a bath with the same temperature, the 
extra energy of the expanding region would have represented 
the conversion to mechanical energy of heat without any 
accompanying work done on the system. However, this is not 
the case with the shock solution. The system is not in ther- 
modynamic equilibrium with any thermal bath or radiation 
field. The temperature is determined only by the balance of a 
heating and cooling rate due to nonthermal processes which 
do not have the same temperature as that assumed for the 
cloud. In the case of the EWS, the pure compression in the 
flowing region of the cloud implies that no extra energy is 
taken in from the heating sources leading to a purely negative 
total energy. The motions of gas in the shock solution are more 
complicated. Initially, static gas is first compressed by the 
passage of the shock. The gas radiates compressional heat, but 
because this gas now has outward motion, it also cools by 
expansion. At some point, the expansion causes the tem- 
perature of the gas to want to dip below the equilibrium tem- 
perature appropriate to the heating and cooling rates. The gas 
then takes in extra energy from the heating sources to remain 
isothermal. Finally, as the parcel of gas is decelerated and sent 
into infall, compressional energy is again radiated by the gas. 
Because there is a net energy taken in by the gas, the energy of 
the flowing region is positive and increasing. 

In the realistic setting of a finite-sized cloud core, collapse 
described by the shock solution would eventually make the 
total energy of the core positive, hence making it unbound. 
This is acceptable behavior in light of the instability of the 
singular isothermal sphere to both collapse and expansion. 
The shock solution simply represents a mode of collapse which 
allows the accretion of some material onto a core but which at 
the same time interacts with an external energy source so as to 
make the flowing parts of the gas unbound. To emphasize this 
point, we demonstrate in § 4 that the singular isothermal 
sphere can be easily sent into complete expansion if collapse is 
not initially forced on the cloud. 

We note that the above energy considerations apply equally 
well to the shock found in the numerical solutions of § 2.2. The 
important thing to consider in that case is that the shock per- 
sists long after the piston has been applied to the gas. The 
energetics of the postshock gas are only related to the actions 
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of the piston in the sense that the piston puts the gas into a 
state where the subsequent evolution is driven by a shock and 
the requirement that the gas be isothermal. The work done by 
the piston is insignificant compared to the subsequent energy 
exchanges of the cloud in trying to remain isothermal. 

We have assumed that the rate at which the cloud is heated 
in the expanding part of the flow is sufficiently rapid so that 
there will only be small departures from isothermality. A 
concern might be that the evolution of the gas in the shock 
solution may require a heating rate which is so high that the 
prescribed processes cannot supply the needed power. This is 
not likely the case since the shock in the above solution is 
weak; the postshock motions are subsonic, and the density 
contrast at the shock is only a factor ~ 1.6. We examine this 
quantitatively. During the early phases of collapse as described 
by the shock solution, the shock is in the high-density region 
near the core. This part of the collapse is not well modeled by 
any isothermal collapse solution since the gas is expected to be 
hotter due to trapped compressional heat. At somewhat later 
times, the shock will be in a more optically thin region where 
grain cooling will be more effective and where isothermal solu- 
tions may apply. Heating in this region can be due to absorp- 
tion of accretion luminosity, cosmic-ray impacts, and grain 
heating. The latter two processes are expected to be most 
important in the outer, low-density regions of the cloud, 
whereas the effects of accretion luminosity are most important 
in the inner regions (although significant heating from accre- 
tion luminosity occurs even when the gas is optically thin; 
Hollenbach et al. 1995; Ceccarelli, Hollenbach, & Tielens 
1995). In Figure 8, we show estimates of heating rates due to 
the aforementioned processes as well as the rate of adiabatic 
cooling due to the expansion of gas immediately behind the 
shock. (We do not include other forms of cooling since these 
are ineffective at temperatures of less than about 10 K.) The 
cooling rate is given by 


A. d (r) = 


5 kTu _ 
2 pm p 


-5 


a 2 up 


(24) 


where all fluid variables are dimensional. We have assumed 



Fig. 8. — The magnitude of the rate of adiabatic cooling of gas expanding 
immediately behind the shock is given by the heavy solid line as a function of 
the position of the shock. The heating rates per unit volume of the same gas 
due to cosmic rays (dotted line) and dust heating (short-dashed line) are also 
plotted as a function of the location of the shock. The long-dashed line indi- 
cates the heating rate due to accretion luminosity. 


that the density distribution obeys p ocr 2 in obtaining the 
rightmost equality. We estimate the heating rates using 


T cr — 4 x 10 28 | — j ergs s 1 cm 3 (25) 

\/""p/ 

for cosmic rays (Black 1987) and 

r dust - 4x 10 _26 (-^-) ergs s 1 cm 3 (26) 

for photoelectric heating from dust (Black 1987). Heating from 
accretion luminosity is estimated by taking the luminosity to 
be 


GM+M 


3 x 10 


"( 


M 


10 5 M 0 yr“ 


1/3 


cm 


(27) 

where the latter relation for R+ is a fit to detailed calculations 
(Adams & Shu 1985; Stahler, Shu, & Taam 1980; Stahler 1983) 
and Af is given by eq. (1) with C = 0.105, the value appropriate 
to the shock solution. The mass of the star is just = Mr, 
where t is the time under consideration. The heating due to 
accretion luminosity for gas at radius r is then computed : 



where k is the absorption opacity of the gas. The dust photo- 
sphere which reprocesses the radiation from the accretion 
shock should have a temperature -500 K (Adams & Shu 
1985) which corresponds to a peak wavelength of —6 pm. The 
opacity here is k — 2.5 cm 2 g" *. Ideally, k should be averaged 
over the emergent spectrum; however, we simply take the 
above value for our estimates. Equation (28) also represents an 
approximation since the emission from the dust photosphere 
will be further reprocessed by the material between the photo- 
sphere and the radius under consideration, r. If r is reasonably 
small, then this should not be too bad an approximation since 
most of the energy fed in at the dust photosphere will even- 
tually make it out to radius r. 

The abscissa in Figure 8 is the location of the shock. The 
heating and cooling rates are computed for the gas imme- 
diately behind the shock because this is the region experiencing 
the greatest expansion. In the inner region (r — 10 15 -10 16 5 
cm) where both cosmic-ray and dust heating are minimized by 
the high optical depth, accretion luminosity provides an ample 
supply of energy to overcome the cooling due to expansion. In 
the outer region (r > 10 16 5 cm) where the amount of accretion 
luminosity absorbed by the gas is minimal, dust and cosmic- 
ray heating adequately replace the energy density lost via 
expansion. 


3.3. Radial Stability of the Shock Solution 

We examine the stability of the shock solution numerically. 
Although less satisfactory than linear stability analysis, essen- 
tial features of the stability of the solution can be discerned in 
this manner. For example, when numerical calculations which 
begin with core collapse are continued beyond the point where 
the EWS is recovered, the subsequent evolution of the gas does 
not deviate substantially from the EWS (e.g., Boss & Black 
1982). This suggests that the EWS is stable to radial pertur- 
bations, a result also determined analytically (e.g., Silk & Suto 
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1988; Ori & Piran 1988). We have also confirmed this result. 
We start a numerical calculation by assuming that the initial 
density and velocity profiles are those of the EWS scaled to 
some early time. The cloud is then evolved for about a dynami- 
cal time with the inner boundary condition that there be free 
fall onto a sink cell (the same boundary condition as adopted 
in § 2). We find no substantial deviation from the EWS. 

The stability of the shock solution can be similarly tested. 
We set the initial density and velocity profiles of a numerical 
calculation to be that of the shock solution at some early time. 
We evolve the cloud with the same free-fall inner boundary 
condition as used previously. The subsequent collapse rapidly 
departs from that of the shock solution and suggests that it is 
unstable. We further test the behavior of the numerical col- 
lapse. Since we know from the analytic solution what the gas 
velocity should be at every location and time, we find that the 
infall speed of the shock solution at the radius where the inner 
boundary condition was imposed is somewhat smaller than the 
corresponding free-fall speed based on the enclosed mass. This 
is because pressure effects, though small, still contribute to the 
dynamics of the flow. In trying to eliminate all possible pertur- 
bations on the flow, we ran another numerical collapse calcu- 
lation imposing the new inner boundary condition that the 
velocity and density be exactly as specified by the analytic 
solution at the given time and location. Even under these con- 
ditions, numerical noise is sufficient to drive the subsequent 
evolution of the gas away from that specified by the shock 
solution. 

The difference in the stability properties of the EWS and the 
shock solution are not surprising. Within the rarefaction wave, 
the EWS possesses a large region of supersonic infall and steep 
velocity gradients which can serve to shear out perturbations. 
This stabilizes the flow (Shu 1977). The shock solution, on the 
other hand, has a large region of subsonic outflow which is 
accompanied by rather limited velocity gradients. Further- 
more, the constant amplitude shock is maintained by a pre- 
carious interplay between a pure outgoing shock and a pure 
outgoing rarefaction wave. Perturbations to this arrangement 
in the form of sound waves can probably do considerable 
damage. A linear stability analysis of the shock solution will be 
presented elsewhere. 

We studied the ways in which the numerical collapse departs 
from the shock solution by perturbing the flow away from the 
analytic solution in various ways. In the first instance, we again 
start the calculation by taking an early-time version of the 
shock solution as giving the initial density and velocity profiles 
for the gas. We also impose the free-fall inner boundary condi- 
tion. However, we artificially increase the initial mass of the 
sink cell by 5%. This has the effect of gravitationally per- 
turbing the entire cloud. The collapse then proceeds according 
to the properties of the class I numerical solutions of § 2.2. 
More specifically, the central mass accretion rate increases 
rapidly from the value specified by the shock solution and 
approaches in time the value corresponding to the EWS solu- 
tion. The amplitude of the shock decreases with time, and the 
ratio of the distances from the shock to the origin and the 
stagnation point to the origin decreases with time (i.e., the 
region of outflowing gas in similarity space gets smaller). Note 
that the same perturbation imposed on a numerical realization 
of the EWS does not result in departures from the analytic 
solution. In the second case, we run a collapse calculation in 
the same manner as above, but we instead artificially decrease 
the initial mass of the sink cell. We first decrease the mass of 


the sink cell so that the free-fall velocity at the radius where the 
inner boundary condition is applied is equal to the velocity as 
specified by the shock solution at that radius. The mass is then 
decreased a further 5% to run the calculation. The subsequent 
collapse evolves according to the class II solutions of § 2.2. The 
central mass accretion rate decreases from the level specified by 
the analytic shock solution, and the shock persists in time with 
roughly constant amplitude. 

The modes of instability described above provide an expla- 
nation for the existence of the two classes of collapse solutions 
found in the piston experiments of § 2.2. After the collapse is 
launched by the piston and imposition of the inner boundary 
condition, the shock propagates out to a distance where the 
initial scale set by the piston motion becomes much smaller 
than the scale of the region over which the gas is flowing. The 
collapse attempts to go toward the shock solution; however, 
because it is unstable, perturbations imposed by the early 
piston motion nudge the collapse toward one of the two modes 
of instability growth. If the piston does not go out too far and 
does not impose initial conditions too far removed from simply 
having an initial pressure imbalance in the core, the solution 
follows the class I type of collapse and evolves toward the 
nearest stable solution, the EWS. On the other hand, with a 
sufficiently large push, the piston launches the collapse into the 
class II type of behavior, with the solution evolving toward 
lower accretion rates. 

The piston experiments show that, given a reasonable 
amount of initial push by the piston, accretion would continue, 
although at low rates, for the duration of the collapse simula- 
tions. It is, however, interesting to consider the limiting case if 
the calculation were to be allowed to run indefinitely. Accre- 
tion will likely have stopped, and we can ask if there is some 
solution which represents pure outflow. That is, are there con- 
ditions under which the singular isothermal sphere could be 
sent into complete expansion? This is considered in the next 
section. 

4. OUTFLOW SOLUTION 

We initiated the numerical calculations of § 2.2 by pushing a 
piston and then imposing a free-fall boundary condition at the 
end of the push. The boundary condition models a large pres- 
sure imbalance in the center of the cloud which essentially 
forces the inner region into collapse. In this section, we con- 
sider the case where a push on the gas by the piston is followed 
by a no-flow boundary condition, i.e., imposing the condition 
that the gas velocity be zero at the inner boundary. 

The results of one of the numerical calculations are shown in 
Figures 9 and 10. The piston again forms a shock which propa- 
gates in time to large radius with constant amplitude. In con- 
trast to the results of § 2.2, the velocity profile is one of pure 
outflow, and the shock has a higher Mach number of roughly 
1.3. The density is characterized by a postshock region which 
flattens out toward small radius, in contrast to the steeply 
peaked density of the initial static isothermal sphere and to the 
density of the shock collapse solution. Identical results to those 
displayed here are obtained assuming any of a series of piston 
speeds and distances (v p and d p ). In particular, the case shown 
in the figures assumes v p = 0.5a and d p — 2 x 10 14 cm. Recall 
that the case from § 2.2 with the same piston parameters, but 
where the free-fall inner boundary condition was imposed, led 
to an increasing mass accretion rate and a diminishing shock 
(Figs, la, 2a, and 3). In the present calculation, the shock does 
not weaken with time. Other cases where d p is even smaller 
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(representing very small perturbations to the flow) also lead to 
constant amplitude shocks. Therefore, for the pure expansion 
case, there is no bifurcation into two classes of solutions. All 
piston parameters considered lead to the same solution as 
shown in Figures 9 and 10. Furthermore, only very modest 
initial perturbations are required to cause the singular isother- 
mal sphere to expand. 

The character of the numerical results again suggests that we 
seek a similarity solution to describe the late time evolution of 
the expanding solution. We use the formulation of § 3.1. The 
shock is again expanding into a static singular isothermal 
sphere. The shock must then obey the jump conditions of 
equation (16). Our procedure is first to assume a location for 
the shock, x 5h . Starting with values of the density and velocity 
given by equation (16), we integrate the fluid equations (8) and 
(9) numerically to small radius. The resulting gas velocity will 
not vanish in general as x -*■ 0. However, by adjusting x sh , this 
condition can be satisfied. We find that such a solution is 
unique, as was the case with the shock collapse solution of 
§ 3.1. The velocity and density profiles are given in Figures 6 
and 7, respectively, by the light solid line. The shock propa- 
gates outward at a constant speed of 1.34a. 

To better compare the analytic result to the numerical solu- 
tions, we scale the similarity solution to the time of the last 
curve in Figures 9 and 10. The similarity solution is given by 
the dots. The numerical results agree exactly with the similarity 
solution, except of course for the shock which occupies several 
zones in the numerical simulation. This indicates both that the 
similarity solution gives a good description of the flow and that 
the solution is probably stable. 

It is easy to show from equations (8) and (9) that 

dec 

- — ► 0 as x -► 0 . (29) 

dx 

The passage of the shock through the singular isothermal 
sphere rearranges the gas into a flat distribution. The centrally 
peaked nature of the initial cloud is completely disrupted as 
gas is sent outward. Figure 7 shows that — 30% of the expand- 
ing region is included in such a “core.” Furthermore, equation 



Fig. 9. — The gas velocities in the pure expansion case where t> = 0.5a and 
d p = 2 x 10 14 cm are plotted for several times. The leftmost curve corresponds 
to the earliest time of 0.1 x 10 12 s. Going toward the right, the corresponding 
times are 0.53, 1.27, 2.37, 3.89, 5.89, 8.44, 11. 6, 15.4, and 19.9 in units of 10 12 s. 
The dots give the analytic solution of Fig. 6 scaled to the time of the last 
numerical calculation shown. 



Fig. 10. — The gas density in the case of pure expansion where v p = 0.5a 
and d p = 2 x 10 14 cm is plotted for the same times as shown in Fig. 9. The dots 
indicate the analytic solution of Fig. 7 scaled to the time of the last numerical 
calculation shown. 

(5) indicates that the central density decreases as the square of 
time. After a time of 10 5 yr, the central density is ^2 x 10 5 
cm -3 , which is comparable to the density of a molecular cloud 
core (e.g., Goldsmith 1987). The expanding solution then not 
only demonstrates that the singular isothermal sphere can be 
easily sent into pure expansion, but that a very modest central 
perturbation and subsequent subsonic motions of the gas can 
lead to a uniform density distribution throughout the region 
originally occupied by the singular isothermal sphere. 

We again emphasize that the shock in the solution of this 
section is not powered by the piston, which stops after a very 
short time. The piston simply puts the gas into a state where 
the subsequent evolution is governed by the dynamics of the 
shock and the assumption that the gas is isothermal. 

5. DISCUSSION AND CONCLUSIONS 

We have considered the collapse of the singular isothermal 
sphere in the case where infall is initiated by some arbitrary 
process which, in addition to causing a central pressure imbal- 
ance, is accompanied by some input of heat or energy into the 
cloud. The rationale for this is to understand the possible 
effects on the collapse of either preheating by accretion lumi- 
nosity or some other means of early energy input into the 
cloud, such as stellar winds. In order to be rather general and 
to understand the relevant physics, we modeled this process 
numerically with a one-dimensional hydrodynamical code 
where we initiated collapse with a central piston which first 
pushes at a constant speed out to some distance which is small 
relative to the size of the cloud. The piston is then stopped and 
a free-fall boundary condition is applied. This inner boundary 
condition acts to initiate collapse in the central regions. We 
assume isothermal conditions throughout. 

Numerical collapse solutions in the piston experiments fall 
into two classes, depending upon how far the piston is allowed 
to push. Both classes are characterized by the presence of an 
initial shock created by the action of the piston, a region of 
outflowing gas just inside of the shock, and inflowing gas in the 
central regions. If the piston is pushed to less than a certain 
distance, the subsequent rate of mass accretion onto the core 
increases with time and approaches the constant value corre- 
sponding to the EWS. The shock decays with time and 
becomes a sonic pulse. The amount of outflowing gas also 
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decreases with time. These were dubbed class I solutions and 
correspond to cases where, because the initial central pertur- 
bation was not very different from the conditions giving rise to 
the EWS, the subsequent collapse tended toward the EWS. In 
the cases where the piston was allowed to push further out, the 
accretion rate instead decreased with time. The initial shock 
persisted with roughly constant amplitude, and the region of 
outflowing gas grew with time. These were dubbed class II 
solutions. These two classes of solutions existed for any of the 
assumed piston speeds which included both subsonic and 
supersonic cases. 

The interpretation of these results is given in terms of the 
instability of a similarity solution (the shock solution) which 
serves to divide the two classes of numerically determined solu- 
tions. The shock solution is characterized by a constant central 
mass accretion rate whose value exactly divides the class I 
accretion rates from the class II accretion rates (Figs. 3 and 5). 
The shock solution also has a constant amplitude shock and a 
constant ratio of the distances from the shock to the origin and 
the stagnation point to the origin. We examine the stability 
properties of the shock solution by starting a numerical col- 
lapse calculation with gas densities and velocities given by the 
shock solution scaled to some early time. The subsequent evo- 
lution of the cloud was found to always deviate from the analy- 
tic solution. (This contrasts with the behavior of the EWS.) 
This indicates that the solution is most likely unstable to radial 
perturbations. The manner of the deviation was such that a 
small perturbation of the flow toward increasing accretion 
rates (artificially increasing the mass of the core), resulted in a 
solution with the exact properties of the class I solutions. A 
perturbation in the other direction (artificially decreasing the 
mass of the core) resulted in evolution akin to that of the class 
II solutions. The interpretation of the numerical solutions is 
that when the shock has reached a distance far removed from 
the scales set by the actions of the piston, the solution tries to 
go to the shock solution. However, because it is unstable, 
initial perturbations due to the piston cause the solution to 
evolve away from the shock solution in either the manner of 
the class I or class II solutions. 

The shock solution and the class II numerical solutions 
eventually lead to the unbinding of any finite region of the 
singular isothermal sphere. The energy for this does not come 
from the actions of any piston. Rather, in the case of the 
numerical solutions, the piston changes the state of the cloud 
so that it absorbs extra energy from whatever agent is helping 
to keep the cloud isothermal. In the case of an actual protostel- 
lar cloud core, this would presumably include dust and cosmic- 
ray heating. Furthermore, this mode of collapse does not 
represent a violation of the instability properties of the singular 
isothermal sphere since it is unstable both to collapse and to 
expansion. The limit of the class II solutions is in fact rep- 
resented by the case of zero mass accretion and pure outflow. 
We also obtained such a solution both numerically and ana- 
lytically. The numerical realization is obtained by again 
pushing a piston out some distance into an initially static sin- 
gular isothermal cloud and imposing the inner boundary con- 
dition that the velocity be zero at small radius. This expanding 
solution is characterized by an outgoing constant amplitude 
shock and a flat central density profile. In approximately a 
dynamical time, the central density falls to values comparable 
to those of interstellar cloud cores. We find that even when the 
piston is pushed a very small distance, a purely expanding 
solution is always obtained. The numerical solution then 


evolves rapidly toward a similarity solution. Because the late 
time evolution of the numerical solutions does not deviate 
from the similarity solution, it is likely that the expanding 
solution is stable. 

We consider the implications of the above results for star 
formation. First, the deposition of energy or pressure into the 
central regions of a collapsing protostellar cloud must occur to 
some extent. For example, if we adopt the standard paradigm 
of low-mass star formation (Shu et al. 1993), as infall is initiated 
in the center of a singular isothermal sphere, a core develops 
rapidly onto which material is accreted. The preheating effect 
of the accretion luminosity will modify the flow in the region of 
gas that is not yet in supersonic infall on a timescale very short 
compared to the dynamical time. This has the effect of giving 
the gas an outward push. Note that the preheating effect does 
not have to be so extensive as to include the more optically 
thin regions of the cloud. If preheating is sufficient to create a 
weak shock in the inner optically thick region, our consider- 
ations show that the shock can propagate to large radius 
simply because the gas is being kept relatively isothermal in the 
thin region. It is also known that protostellar winds occur very 
early on in the formation of the protostar (e.g., see Bally & 
Lane 1991). Although these winds are thought to be initially 
highly collimated, they nevertheless can impart significant 
outward momentum to some of the gas. We have modeled all 
of these processes by a piston pushing in the center of a singu- 
lar isothermal sphere. Although heuristic, this approach does 
give an indication of the nature of the subsequent collapse. 
Specifically, there is the possibility of accretion rates that are 
lower than the canonical value given by the EWS. Depending 
upon the degree of the initial perturbation, accretion rates 
could be a factor of 10 lower than the EWS value and evolve 
toward lower values, or they could be initially rather low and 
rise up to the EWS value. The present work does not predict 
the exact rate since the extent of the initial perturbation is not 
known. 

The prediction of reduced accretion rates may be related to 
the “luminosity problem" of protostars. The observed lumi- 
nosity function of embedded sources in several star-forming 
regions implies that the central accretion rates are generally 
~ 10 times lower than predicted by the standard paradigm 
(e.g., Kenyon et al. 1990; Kenyon et al. 1994; Greene et al. 
1994). This discrepancy may be due to rotational effects that 
are not considered in the spherically symmetric approximation 
of the EWS (Kenyon et al. 1994). For example, most collapsing 
material could fall onto a disk and then be accreted onto the 
star on a timescale set by disk physics. The accretion from a 
disk could be, for the most part, smaller than that of spher- 
ically symmetric accretion. Since studies show that the lifetime 
of the accretion phase is ~ 10 5 yr, accretion from the disk must 
occur in a series of bursts separated by a phase of low accretion 
in order that the average accretion rate be high enough to form 
a solar-mass star in the given time. Our study indicates that 
independent of rotational considerations, which must be 
important to some degree, spherically symmetric collapse itself 
gives rise to time variable accretion rates which could be rather 
low, especially at early times. The objection is then raised that 
in ~ 10 5 yr a consistently low accretion rate leads to a star that 
is well below a solar mass. There are several possible ways out 
of this. A time- variable accretion rate, such as those of the class 
I solutions, would imply that the low-luminosity sources are 
simply young and are experiencing a period of rather low acc- 
retion. At later times, the accretion rate increases so that a 
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solar-mass star can still be formed in approximately the 
allowed time. On the other hand, if indeed the accretion rate 
remains consistently low, in ~ 10 5 yr the shock of our solutions 
would have traveled to the outer regions of the cloud core. The 
nature of the collapse after this depends in part on what 
happens to the shock on reaching the outer boundary, 
although the presence of significant magnetic support for the 
cloud in these regions will lessen the effects of the shock. Since 
this is unknown, what happens to the accretion rate after the 
first ~ 10 5 yr is unknown also. 

The application of our study to protostellar winds is more 
heuristic. Early stellar winds are thought to be highly colli- 
mated and rather supersonic. Our solutions, however, are 
spherically symmetric, and the shock is rather weak. Therefore, 
our considerations will not provide an exact description of the 
effects of the wind on the collapsing cloud. However, we 
showed that the singular isothermal sphere is, qualitatively 
speaking, rather loosely bound; small initial perturbations 
toward expansion can easily lead to the general expansion of 
the cloud. This indicates that the termination of accretion can 
be easily effected by a stellar wind. Conversely, the star forma- 
tion scenario might be that protostellar collapse first proceeds 
according to one of our solutions. The passage of the shock 
unbinds much of the material in the outer regions of the core. 
The onset of an observable protostellar wind may then occur 
as a result of the cloud no longer presenting much of a barrier 
against which the wind must fight to get out of the core. 

The results of § 4 also offer the intriguing possibility that if 
some process exists which applies a small push in the central 
regions of a centrally peaked cloud core without an accom- 
panying pressure imbalance, the core can be easily dispersed in 
a dynamical time. It is not clear, however, whether such a 
process exists in nature. 

In considering the relevance of the current study to star 
formation, we must be mindful of the limitations of our 
assumptions and methods. For example, the assumption that 
protostellar collapse occurs in an isothermal fashion has pre- 
viously only been made in connection with a purely collapsing 
cloud. Under this latter condition, compressional heat gener- 
ated during collapse is efficiently radiated by various mecha- 
nisms, including dust emission. However, our solutions include 
an expanding outer envelope which is cooled by expansion. We 
should therefore check that sufficient heating exists from 
various processes (cosmic rays, dust heating, and heating from 
accretion luminosity) to keep the gas isothermal. We have 
demonstrated that this is indeed the case by estimating the 
various heating rates and comparing them to the rate of adia- 
batic cooling of expanding gas (see § 3.2). 

Another assumption of the current work is that the cloud 
into which the shock propagates is exactly stationary. If the 
currently held paradigm for low-mass star formation does 
indeed provide a framework for collapse (Shu et al. 1993), there 
could be a subsonic inward drift of the gas at large radius as 
the magnetic field diffuses out of the collapsing cloud. What is 
the behavior of the shock solution under these conditions or in 
the more extreme case where the inward drift is supersonic? It 
is obvious, since the shock in our solution has a Mach number 


of only 1.26, that any inflow upstream of the shock of greater 
than the indicated speed will sweep the shock into the core. 
The solutions investigated in this paper would not apply to 
these flows. In the case where the inflow ahead of the shock is 
subsonic, the effects of our solution are somewhat reduced. For 
example, we have considered the collapse of isothermal spheres 
which have a density that is everywhere greater by a uniform 
factor than that of the singular isothermal sphere. Because 
these spheres are too dense to be in equilibrium, infall is initi- 
ated at all radii at the beginning of collapse. The case where 
infall is not accompanied by any central energy input is 
described by the MSWCP solutions discussed in § 3.1. When a 
central piston is applied, we find that resulting solutions are 
similar to those in which a shock propagates into static gas; 
however, there is a greater tendency to follow the class I behav- 
ior (see § 2.2). Specifically, the accretion rate, which always 
begins at rather low values (M ~ few x 1CT 7 Af 0 yr ! ) shows 
a greater tendency to evolve toward the higher value associ- 
ated with the corresponding MSWCP solution. The low accre- 
tion rates characterizing the early parts of the collapse of some 
of the class I solutions last for shorter times. Similarly, the 
magnitude of the shock tends to decrease more rapidly as does 
the size of the region of outflowing gas. However, even when 
gas ahead of the shock is infalling subsonically, significant 
reductions in the central accretion rate relative to the values 
associated with no energetic input are possible. 

Current observational constraints on the nature of infall are 
inadequate for differentiating between our solutions and other 
models. Although molecular line width studies have ruled out 
outflow in various cloud cores (e.g., Zhou 1992; Evans et al. 
1992), the nature of the infall has not been especially well con- 
strained. Although collapse according to the Larson- Penston 
model (Larson 1969; Penston 1969a) is likely ruled out, the 
nature of the infall predicted by the shock solution is not very 
different from that of the EWS, in that free-fall will prevail in 
the central regions. Additionally, the outflows in our solutions 
are rather subsonic, so they are not detectable using line pro- 
files. The constraints on the density profile are similarly ambig- 
uous enough to allow a variety of models. Studies of spectral 
energy distributions indicate the density can be consistent with 
p ~ 5-2.0) p ro fii es (Butner et al. 1990; Butner et al. 1991). 

Because the EWS and our solutions both derive from collaps- 
ing a singular isothermal sphere, comparable density profiles 
exist throughout much of the cloud. 
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